Fcc-bcc transition for Yukawa interactions determined by applied strain deformation 
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Calculations of the work required to transform between bcc and fee phases yield a high-precision 
bcc-fcc transition line for monodisperse point Yukawa (screened-Couloumb) systems. Our results 
agree qualitatively but not quantitatively with recently published simulations and phenomenological 
criteria for the bcc-fcc transition. In particular, the bcc-fcc-fluid triple point lies at a higher inverse 
screening length than previously reported. 
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I. INTRODUCTION 

The screened- Coulomb or Yukawa pair potential, 
U(r) = ^e~^^ /r, has been the focus of great theoreti- 
cal interest for two reasons. One is that it describes a 
wide range of interactions, changing continuously from a 
pure Coulomb potential to an effective hard-sphere po- 
tential as the inverse screening length k increases. The 
second is that it provides an approximate description 
of the effective interactions between large ions that are 
screened by more mobile counterions. In this context 
it has been used to describe the interactions between 
ions surrounded by electrons in metals Q .dust grains 
surrounded by electrons in dusty plasmas [1, 0, 0| , and 
macroions surrounded by counterions in charge-stabilized 
colloidal suspensions 0, IE IE S Hi • 

The phase diagram of systems of particles interact- 
ing with a Yukawa potential has been studied with 
both analytic 1^ liol [ill IT^ flsL 113 and numerical 
HI EE E 111 COal^niilllllll techniques 
and compared to experiments on dusty plasmas ^] 
and colloidal suspensions IH IIE • The high temper- 
ature phase is a fluid. There is no liquid-gas transition 
because the interactions are purely repulsive. The sta- 
ble crystalline phase at zero temperature changes from 
bcc to fee as K increases. The higher entropy of the bcc 
phase leads to a greater range of stability as temperature 
increases until the melting line is reached. Previous re- 
sults for the fcc-bcc transition line [11II20L I21II2I I2I |23 | 
vary substantially and the most recent detailed calcu- 
lation ^El quotes an uncertainty of about 10% roughly 
halfway between the zero-temperature transition point 
and the triple point. 

In this paper we use a different approach to obtain 
the bcc-fcc phase boundary with an uncertainty of only 
about 1%. Bounds on the free energy difference between 
the two phases are obtained by calculating the work done 
during a continuous deformation between them. The ef- 
fect of deformation rate, truncation of the potential, and 
system size and geometry are all analyzed to determine 
systematic errors. The resulting bcc-fcc transition line is 
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in qualitative agreement with recent simulation results, 
and quantitative differences are comparable to the larger 
error bars quoted by previous studies. We estimate the 
location of the bcc-fcc-fluid tri ple point using previously 
published melting-line results [TtI l23j| , and find that it 
lies at higher inverse screening lengths than previously 
reported. 

The role of anharmonicity in stabilizing the fee phase 
is analyzed in detail. While anharmonicity increases the 
energy of the fee phase relative to that of the bcc, there is 
an even larger increase in the relative entropic contribu- 
tion to the free energy that increases the range of stability 
of the fee phase. This appears to reflect an increase in 
the frequency of the long wavelength shear modes that 
dominate the bcc entropy in the harmonic approximation 

Our results are also compared to phenomenological cri- 
teria proposed by Vaulina et. al. [iE|- These authors pre- 
dict a transition at a critical value of the mean-squared 
displacement about lattice sites, and calculate the dis- 
placement using a simple Einstein-like model. We find 
that the actual displacement from MD simulations on 
our transition line is in reasonable agreement with their 
phenomenological criterion, but substantially larger than 
predicted by their Einstein model. 

The details of our calculations are presented in the fol- 
lowing section. Section ITTll provides a detailed analysis of 
systematic errors and presents our results for the phase 
boundary. In Section llVI we compare our results to pre- 
vious transition lines, and SectionlVlprovides a summary 
and conclusions. 



II. METHOD 

A. Free energy difference calculations 

NVT ensembles are most natural for the study of 
Yukawa systems for two reasons. First, since the Yukawa 
potential is purely repulsive, the macroions in an exper- 
iment will expand to fill the container. Second, the in- 
verse screening length k is density dependent in charged 
colloidal suspensions and dusty plasmas IsJ. This den- 
sity dependence is system-specific, and affects the pres- 
sures and bulk moduli. Thus any calculation of coex- 
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/•bcc 

Wft = V a- de, 

J fee 



(2) 



where a and e are the stress and true strain tensors. Wbf 
and {—Wfb) are upper and lower bounds on AF. For 
these bounds to be narrow, the intermediate configura- 
tions of our systems must remain statisticahy represen- 
tative of the ^-dependent equihbrium distributions as C 
is varied [30|. In particular, the stress tensor ct(C) must 
remain near its equilibrium value. 



B. Potential parameters 




FIG. 1: The Bain transformation. Two cells of a bcc lattice 
are shown with lattice directions. The body-center atoms are 
shown in black. When the x, y, and z directions are scaled by 
2-1/6^ 2^/^), the crystal is transformed into an fee lat- 
tice of the same density. The atoms connected by the dotted 
lines become two (100) faces of an fee unit cell. 



istence regions will be non-universal. For this reason 
we focus on finding the Helmholtz free energy difference 
AF = Ffcc ~ Fbcc at fixed volume. A brief discussion of 
coexistence is given in Section FlI Bl 

Postma, Reinhardt, and others [sollsills^ have shown 
that the free energy difference between two phases of a 
system may be calculated in numerical simulations by 
evaluating the external work done on the system along 
a thermodynamic path connecting the phases. From el- 
ementary thermodynamics, the mechanical work Wab 
done on a system on an isothermal path from state A to 
state B gives an upper bound on the change AFab = 
Fb — Fa in the system's Helmholtz free energy jS^ • The 
work Wba done on the system during the reverse process 
B ^ A is an upper bound on AFba, and hence {—Wba) 
is a lower bound on AFab- 

Bounds on AF for Yukawa systems can be obtained us- 
ing a continuous constant-volume Bain deformation path 
(Fig. connecting the bcc and fee lattices. An initially 
bcc lattice deformed such that its three cubic symmetry 
directions are scaled respectively by X) is 

transformed into an fee lattice of the same density as C 
varies continuously from 1 to 2^/^ calculate the 

work done along this path in the forward and reverse di- 
rections using strain-controlled molecular dynamics sim- 
ulations [s^. 

Assuming that the systems traverse these paths ho- 
mogeneously, we can calculate the work done from the 
global stresses and strains. We define 



The phase behavior of Yukawa systems is most con- 
veniently expressed in terms of dimensionless screening 
and temperature parameters. One natural, phase- inde- 
pendent length scale is a = n~^/^, where n = N/V is the 
macroion number density. The Yukawa potential may 
then be expressed as 



U{r) 



a r/a ' 



(3) 



ffcc 

Wbf^V a -de, 

J bcc 



(1) 



where X = na is the dimensionless screening parame- 
ter. The limits A ^ and A — > cxd correspond to the 
exhaustively studied one-component plasma and hard- 
sphere systems. 

A natural time scale is provided by te, the period of an 
Einstein oscillator in a crystal. The Einstein periods for 
the fee and bcc phases change by an order of magnitude 
over the range of A studied here (3 < A < 8), yet differ 
from each other by less than 1.2% at any given A within 
this range. To obtain consistent results across a wide 
range of screening lengths, we normalize all time scales 
in this study to te{X), using the fee values given in Ref. 

M 

A natural energy scale is given by the Einstein phonon 
energies muP^a? , where m is the macroion mass and 
loe — is the Einstein frequency. Following Kre- 

mer, Robbins, and Grest |20[, we define the dimensionless 
temperature 

(4) 

using the fee phonon energies, and plot our phase dia- 
gram in (A, T) space. A dimensionless inverse temper- 
ature F = {^/a)/kBT called the coupHng parameter is 
used in many studies of dusty plasmas. The advantage 
of using T rather than F in Yukawa phase diagrams is 
that the transition lines are approximately linear in A. 

The bcc and fee phases coexist in equilibrium over a 
par t of the phase diagram. Following previous authors 
|20l I21L I23 . 12^ , we define the bcc- fee transition line as 
the curve ftrans{\) on which AF = AF(A,f ) ^ 0. This 
transition line will certainly lie within the coexistence re- 
gion, regardless of the thermodynamic state dependence 
of K and We find Ttrans{^) by calculating AF at many 
points {\i,Ti) on the phase diagram. 



3 



C. MD simulation details 

We simulate NVT ensembles of identical particles us- 
ing a velocity- Verlet [S^ algorithm to integrate the par- 
ticle trajectories. The temperature is maintained with a 
Langevin thermostat |37j |. Periodic boundary conditions 
are used to maintain the density. The equations of mo- 
tion for the position qi and peculiar momentum pi of the 
ith particle are 

Pi ^ Fi~ epi +ffi- Pi/TLang, 

where e is the true strain rate tensor, Fi is the force due 
to Yukawa interactions, ffi is a random noise term, and 
TLang IS the characteristic relaxation time of the ther- 
mostat. We use a timestep St = .OIte to insure proper 
integration of Eqs. |(SJ) and set r^ang — IOt^. Changing 
6t and TLang by a factor of two in either direction had no 
effect on the phase diagram. 

For numerical efficiency we truncate interactions at a 
cutoff radius Tc- Due to the presence of long range order 
in Yukawa crystals, care must be taken in choosing this 
cutoff radius. We present the details of our determination 
of rc(A) in Section UTTbI 

In most of our simulations, we impose the bcc^fcc 
Bain transformation as follows. We start with a lattice 
of 3456 particles (12'^ bcc unit cells) in a cubic simu- 
lation cell with edges of length = Ly = Lz = Lq 
aligned with the < 100 > directions of the lattice. The 
system is equilibrated for 200 Einstein periods. We then 
fix ^ = Lz/Lq for a time At sufficient to reach the fee 
structure: At = (23 — The other cell edges Lx and 

Ly are varied to maintain constant volume and tetrago- 
nality (L^ = Ly = \/ Li^/ L^) ■ The true strain rate tensor 
e is then given by 

ezz = Lz/Lz^ t/C,, 

^xy — ^xz ^yz 0- 

We compute the diagonal elements [Px, Py, Pz) of the 
pressure tensor using standard methods 38]. Equation 
then takes on the more physically familiar form 

Wbf ~ {PxLyLzLx + PyLxLzLy + PzLxLyLz^ dt . 

Jo 

(7) 

After the system has reached the fee structure, the de- 
formation process is reversed by changing the sign of 
As the system returns to bcc, Wfb is calculated using the 
analogue of Eq. 

To minimize uncertainties in Titans (A), C must be small 
enough for the system to remain near equilibrium. One 
requirement is that the strain-rate components {eqi) of 
the velocities must be small compared to the thermal ve- 
locity. The Bain transformation time At (which is pro- 
portional to C~^) must also be large compared to TLang to 



allow the thermostat to transfer heat to or away from the 
system as necessary to maintain constant temperature. 
Since T^ang sets the time over which the system samples 
the canonical ensemble, the thermodynamic sampling im- 
proves as At /TLang iucrcascs. 

The precision of the calculated transition line depends 
on the difference {AFmax - AFmin) = Wbf + Wfb = 
Wcycie between the bounds on AF. These bounds con- 
verge to each other in the reversible thermodynamic 
(zero strain rate) limit. In this limit, the average work 

< Wcycie > done on the system over a full deformation 
cycle (bcc— >fcc^bcc or vice versa) should vanish. In sim- 
ulations at finite strain rate, however, there is a positive 
systematic error in Wcycie due to energy dissipation |39j . 
This can be physically interpreted as arising from vis- 
cosity. Each applied strain increment takes the system 
slightly out of equilibrium. When Q is small, one expects 
the stresses to deviate from their equilibrium values by 
an amount avisc ^ C I^O]. Sources of viscous dissipa- 
tion include the intrinsic viscosity and the drag forces 
—pi /TLang On the particlcs applied by the Langevin ther- 
mostat. The viscous dissipation rate is given by Umsc ■ e, 
so one expects the dissipated power to be proportional 
to C/^ . Since the total simulation time scales as C~^, 
the total dissipated energy, and hence the deviation of 

< Wcycie > from zero, should be linearly proportional to 
C. We present the C-dependence of our results in Section 

inrsi 



III. RESULTS 
A. Strain rate dependence 

At temperatures near the transition line, calculations 
of the geometrical structure factor and pair correlation 
function verify that our systems traverse the Bain trans- 
formations homogeneously. This homogeneity allows us 
to use Eqs. H1I2|I for calculating Wbf and Wfb and leads 
to tight bounds on AF. 




-3 t 



FIG. 2: Simulation results for A = 5, C = IO'^/te. Solid 
triangles are values of Wbf, solid squares are values of Wfb- 
The dashed and solid lines are linear fits to the results. Tbf 
and Tfb are the intersections of these lines with W = 0. 
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FIG. 3: Strain- rate dependence of 5hyst- Solid circles indicate 
A = 4 results. Empty squares indicate A = 5 results. The 
solid and dashed lines are linear fits to the data. The error 
bars indicate statistical uncertainties. (Color online) 



We obtained results similar to those shown in Figure 
El over the entire investigated range of A and for three 
different values of Q. Near the transition line, Wbf{\, T) 
and Wfb{\T) vary linearly with T and have nearly op- 
posite (A-dependent) slopes. The scatter about hnear 
fits to W),}{\, T) and W/f,(A, T) is consistent with fluctu- 
ations in Whf and Wfb at fixed (A,T). The intersections 
of these fits with W — Q give two estimates, Tf,/ and T/;,, 
for Ttrans- Thcse are obtained using data at ten evenly 
spaced T within about 5% of the transition line. 

For a given system, Tjb and Tbf provide upper and 
lower bounds on Ttrans since Wbf > AF > ~Wfb and 
(dAF/dT) > 0. We define the fractional uncertainty 
due to dissipative hysteresis as 



Shyst - [Tfb - Tbf)/ {Tfb + Tbf). 



(8) 



The conditions 6hyst ~ C and < Wcycie >'~-" C are equiva- 
lent due to the linear dependence of Wb / and Wfb on T. 
Figure 13 shows the ('-dependence of Shyst for A = 4 and 
A = 5. The results are consistent with our hypothesis 
that the energy dissipated is linear in (. 

We identify Ttrans — {Tbf+Tfb)/2 as the best estimate 
for the transition temperature for a given system size, 
system geometry, and potential cutoff radius. Table HI 
shows that the fractional variation of Ttrans with ( is 
much smaller than Shyst 113 ■ 



TABLE I: 


Dependence of Ttrans (A) 


on dimensionless strain 


rate. 






(te 


10^ftra„4A = 4) 


W'^Ttransi^ ~ 5) 


5-10-5 


1.637 ± .003 


2.363 ± .002 


1 • 10-'' 


1.634 ± .004 


2.365 ± .004 


2 ■ 10-" 


1.633 ± .005 


2.366 ± .005 



In the following, we present results for \(\ = IQ-^/rs. 
Based on Tabled for this value of \^\ the random and 
finite strain rate uncertainties in Ttrans are comparable, 
both about 0.2%. The combined error is estimated to be 
less than 0.4%. The uncertainties given in subsequent ta- 
bles include only statistical uncertainties from the linear 
fits used to calculate Tbf and Tfb- 

B. Potential cutoff dependence 

We estimate the errors introduced by truncating the 
force at Tc by calculating the error 5E{rc) in the potential 
energy difference. If the error in AF is of the same order, 
then the fractional error in the transition temperature is 



_ 5T{rc) 1 dT 



TdAF 



5E{rc)- 



(9) 



Here (dT/dAF) is known near the transition line from 
the work calculations. 

The cutoff-induced error in the potential energy dif- 
ference can be written in terms of the pair correlation 
functions gfccir) and gbccif) of the fee and bcc crystals. 
For N particles 



N 



U{r)igfUr)-gbcc{r))ATTr^dr/a'' (10) 



If = max{\gfcc{r) 



\5E{r,)\ 
7V($/a) 



U{r)2nr dr 



(r)|; r > rd, then 



A2 



(11) 

One expects f2 to be of order 1 at finite temperature. 

For A = 3, 4, and 5, we estimated SE{rc) at T ~ Ttrans 
by calculating the pair correlation functions in large sys- 
tems using large cutoff radii and long integration times. 
Due to the exponential falloff of U(r) and finite temper- 
ature smoothing of g(r), the infinite upper bound in Eq. 
I10|l can be replaced by a finite value r; without introduc- 
ing significant errors. We found Ar; = 30a to be suffi- 
ciently large. 

Figure 0] shows our estimate of \SE{rc)/ {N^/a)\ from 
Eq. (|10|l and the value of \6E/ {N^/a)\ corresponding 
to Scut = 0.01 for A = 3, the longest-range potential 
considered. We found that decreases from 1.4 to .62 as 

increases from 3.5a to 6.5a. The actual error is always 
smaller than the bound given by because gfcc — gbcc 
oscillates in sign. The envelope shown corresponds to 
0.-1/3. Because of the sharp variation of dE{rc), we 
use the envelope to estimate Scut- 

To test Eq. we calculated Ttrans as a function of Vc 
for A = 3. Results are shown in Tabled] The fractional 
changes in Ttrans from Tc — 3.5a and rc = 4.667a to Tc = 
6.667a are 19% and 0.9%, respectively. Both changes 
are about one-fifth of the estimates for Scut from Eqs. 
(|9I10|) . No statistically significant changes are expected 
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FIG. 4: Determination of for A = 3. The heavy curve 
is our estimate of SE{rc) obtained from simulations at T = 
8.84 ■ 10""^, with rj = 10a. The dashed line is the analytic 
upper bound on 5E{rc) from Ea. pijl with Q. = 1/3. The 
horizontal line is the value of \5E\ in Eq.Q corresponding to 
\5,ut\ =0.01. 



or observed for Tc > 5.833a. We conclude that errors 
estimated from the envelopes of curves like Fig. 01 give a 
conservative estimate of cutoff errors. 



TABLE II; Ttrans vs. for A 
and the bound in Eq. 1111 



3. \Scut\ is given by Eq. El 



3.5a 
4.667a 
5.833a 
6.667a 



\Scut 

1.04 
0.041 
0.0015 
0.00014 



1 Tt rans 

1.048 ± .003 
0.871 ± .002 
0.880 ± .002 
0.879 ± .003 



ically lower, but the effect was small. From theoretical 
considerations one expects the leading finite size correc- 
tions to Ai^/iV to be proportional to l/N [ij. This 
should produce a corresponding error in Ttrans- As shown 
in Table IfffI the changes in Ttrans from N — 432 to 
N — 3456 were all about 1%. The changes in Ttrans from 
N — 3456 to = oo for this system geometry should be 
about 8 times smaller. 



TABLE III: Dependence of Tt, 



on N. 



10=*Tt,„„(iV = 432) KfTtransiN = 3456) 

2.350 ± .009 2.362 ± .004 

2.985 ± .011 3.021 ± .003 

3.562 ± .009 3.593 ± .005 

4.064 ± .008 4.087 ± .005 



A 
5 
6 
7 



Another test indicates that finite size effects are larger 
than the above estimate. The geometry was changed so 
that the fee state has equal cell edges and the bcc state 
has Lx = Ly = \/2Lz. These simulations contained 10'^ 
fee unit cells (4000 particles) with the < 100 > directions 
parallel to the simulation cell edges. After Bain transfor- 
mation, the bcc state has two < 110 > directions parallel 
to the simulation cell edges. As shown in Table HVl the 
values of Ttrans obtained for both A = 4 and A = 7 were 
0.6% lower than those obtained with the standard sys- 
tem geometry |53j . Other simulations verified that this 
was due solely to the change in boundary conditions. We 
conclude that our dominant source of uncertainty is finite 
size and is less than 1%. 



To ensure that the fractional systematic errors were 
no larger than our random and rate errors, we chose Vc 
slightly above the values corresponding to \Scut \ — 0.002. 
For A = 3, 4, and 5, we used cutoff radii of 5.833a, 4.375a, 
and 3.5a in the simulations used to determine Ttrans- 
Smaller Tc can be used at higher A both because the 
interactions weaken and Ttrans/Tmeit increases, leading 
to a smaller $7. For A > 5 we fixed the cutoff radius at 
Tc — 3.5a. 



C. System size and geometry dependence 

To examine finite size effects we also considered a 432- 
particle system (initial state 6^ bcc unit cells). Because 
the corresponding fee state has transverse length 6.73a, 
the minimum image convention requires Tc < 3.367a, and 
we used Tc — 3.3a. To separate out rc-dependence from 
system size dependence, we also recalculated the transi- 
tion line for N = 3456 for 5 < A < 8 with Tc = 3.3a. 

Table Hill shows a comparison of our calculated transi- 
tion temperatures. The N = 432 values were systemat- 



TABLE IV: Dependence of Ttrans on system geometry. 



A 


10^rt™,(7V = 3456) 


10^rt,a,„(iV = 4000) 


4 


1.634 ± .004 


1.625 ± .004 


7 


3.592 ± .004 


3.570 ± .004 



We attribute the observed sensitivity to geometry to 
the change in allowed low frequency modes. These modes 
play a disproportionate role in determining the entropy in 
lattice dynamics calculations (421 and drive the fcc^bcc 
transition with increasing temperature [2]| . Since the 
shear velocity is highly anisotropic in the bcc phase, 
changing the boundaries affects the sampling of these low 
frequency modes and thus Ai^. 



D. Transition line 

TableOshows our calculated Ttrans{^) with statistical 
uncertainties. As described above, the combined system- 
atic errors due to finite strain rate, system size, and po- 
tential cutoff are estimated to be less than 1%. Results 
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FIG. 5: Phase diagram of Yukawa systems. The heavy hne is 
our cubic polynomial fit for the bcc-fcc transition line. The 
light solid line is the fcc-bcc transition line from Ref. . The 
low and high dotted lines are respectively the lattice dynamics 
and molecular dynamics bcc-fcc transition lines from Ref. [2H . 
The triangles are a bcc-fcc coexistence point and triple point 
from Ref. [2^ . and the solid square is a bcc-fcc transition 
point from Ref. [2^ . The dashed and dash-dotted lines are 
the melting lines from Ref. [l^ and Ref. |2^. The empty 
squares denote our estimates for the bcc- fee- fluid triple point. 



pressure is larger than the fee pressure by only about 
0.04% for A = 4, and by 0.65% for A = 7. This results 
in a higher density in the fee phase at coexistence, but 
only by about 0.015% at A = 4 and 0.2% at A = 7. The 
corresponding changes in A and T are much smaller than 
the uncertainties in our calculated transition line. The 
coexistence region in experimental systems may be much 
larger due to variations in k and $ with density (26ii . . 
As noted above, these variations are system specific and 
a more complete treatment is beyond the scope of this 
paper. 



E. Anharmonic effects 

The dotted line in Figure [S] shows lattice dynamics 
results for the bcc-fcc transition line ^21] ■ In this approx- 
imation the energy and entropy diff'erences, A_Bld and 
AS'iD, are independent of T. The fcc-bcc transition line 
is given by T^c = /^^Eld / /!^Sld- The resulting curve 
lies below Ttrans, indicating that the fee phase is stabi- 
lized by anharmonic effects |^ |2^ . This implies that 
the anharmonic component of the free energy difference. 



of a cubic polynomial fit to the data are also given: 

10''t;{.;*„^(A) = 6.46678(A - Aq) + 0.43001(A - Ao)^ 
-0.06806(A- Ao)^ 

(12) 

where Aq = 1.718 is the zero-temperature transition point 
obtained from lattice statics calculations Lower- 
order polynomials fail to adequately fit the data within 
our uncertainties. 

TABLE V: Calculated and fit values of Ttrans- Only statisti- 
cal uncertainties are quoted. 



A IQ^'Ttrans ^^'Tl^i 

3 0.880 ± .002 0.885 

4 1.634 ± .004 1.619 

5 2.365 ± .004 2.345 

6 3.017 ±.004 3.023 

7 3.592 ± .004 3.613 

8 4.085 ± .004 4.072 



Figure shows the polynomial fit and two previously 
published solid- fluid coexistence lines [l^ |2^ . The inter- 
sections of these lines give estimated values of the bcc- 
fcc-fluid triple point. Using results from Ref. ^3 ^^'^ 
(Xtp = 7.45, ftp = 0.00384). Those from Ref. H yield 
(Atp = 7.84, ftp = 0.00401). 

If we assume that the parameters k and $ are density- 
independent, we can calculate the width of the bcc-fcc 
coexistence region from the pressures and bulk moduli 
of the two phases on the line where AF = 0. The bcc 



AFan=AEan-TASan^AF-AFLD, (13) 

is negative on the transition line. The relative signs and 
magnitudes of AEan and ASan may be calculated by 
comparing our accurate measurements of free and total 
energy differences with the lattice-dynamics results. 

Table IVII shows results for anharmonic contributions 
to the free and total energy differences on the fit transi- 
tion line IH^. The values of AFan are known from the 
work calculations, while the values of AEan were found 
from separate equilibrium simulations. The anharmonic 
corrections to the total energy favor the bcc phase for all 
A, i.e., Efcc — Ebcc on the transition line has increased 
relative to its zero-temperature value. The anharmonic 
contributions to the free energy difference, however, are 
larger in magnitude and opposite in sign, implying that 
anharmonic entropic contributions to AF favor the fee 
phase at all A and overwhelm energetic contributions. 

TABLE VL Anharmonic free and total energy differences 
evaluated at f//„*„^(A). 



A W^AFan/NkBT W^AEan/NkBT 

3 -0.58 ±0.1 0.42 ±0.33 

4 -1.07 ±.12 1.0 ±.4 

5 -1.92 ±.15 1.2 ±.4 

6 -3.03 ±.17 2.4 ±.4 

7 -4. 33 ±.21 2. 2 ±.4 



In lattice calculations the larger entropy of the bcc 
phase comes mainly from the lower frequency of its shear 
modes. Some of these modes have negative energy for 
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A > 7.67, causing the bcc pha se to become hnearly un- 
stable at low temperatures |21| . It is interesting that the 
onset of this low temperature instability is close to Xtp- 
However, we have performed runs near the melting line 
for A as large as 10 and find that the bcc phase remains 
metastable. This implies that anharmonic effects have 
increased the frequency of long wavelength shear modes, 
providing an explanation for the decreased entropy ad- 
vantage of the bcc phase. 



IV. COMPARISON TO PREVIOUS RESULTS 

Miller and Reinhardt were the first authors to use Bain 
deformation paths to obtain bounds on AF for Yukawa 
systems . They calculated the work by integrating the 
change in the Hamiltonian rather than from the stresses 
and strains. The large discrepancy between their A = 7 
transition temperature and our result is likely due to their 
extremely small system size {N = 108), which was just 
used to illustrate their method. 

The earliest MD calculations poL l2l| of the transition 
line also deviate substantially from ours, particularly at 
large A. The line shown in Figure [3 is a fit between 
points where the fee and bcc phases were found to be 
stable. The gap between points was about 20% and the 
final shape was strongly influenced by a bcc-stable point 
above the melting line. Other points where the bcc phase 
was stable lie close to our Ttrans but are shifted up due 
to the smaller used. 

Our transition line is in qualitative agreement with 
more recent MD and Monte Carlo results [23. l23l| . 
Dupont et. al. calculated a fcc-bcc coexistence point 
and the fcc-bcc-fluid triple point using small systems 
{N ~ 250). Ahhough their triple point {Xtp = 6.75, ftp = 
0.0034) lies wefl below ours {Xtp = 7.7 ± 0.3, ftp ^ 
0.0039 ±0.0001), it lies only about 2% below our fcc-bcc 
transition line, and well below recently published melting 
lines mill. 

Hamaguchi, Farouki, and Dubin also obtain a lower 
triple point (Atp ~ 6.90, ftp — 0.0038) because their 
bcc-fcc transition temperatures are systematically (6- 
10%) higher than ours [2|. One possible explanation 
is that their equilibration times were too short. They 
used the A-independent time unit r = , where 

Ljp — \J ^-nn^ jm is the plasma frequency. Starting with 
perfect bcc and fee lattices as their initial conditions, they 
equilibrated their systems for a maximum of 300r before 
beginning their free energy measurements. This corre- 
sponds to about TJte for A = 3 and only \te for A = 8. 
Since the latter is only about four times the velocity au- 
tocorrelation time, and comparable to the time for sound 
to propagate across their simulation cells, it is doubtful 
that their systems had equilibrated sufficiently at high A. 
Too short an equilibration time could cause overestima- 
tion of the stability of the phase with lower entropy, the 
fee phase, which is consistent with their findings. 
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FIG. 6: Comparison of transition line to analytic estimates. 
The heavy line is our t/^^^^{X). The dashed and dotted lines 
are the analytic estimates for the bcc— >fcc and fcc-bcc tran- 
sitions from Ref. flo| . 

Vaulina and colleagues have proposed phenomenolog- 
ical criteria for the bcc-fcc transition 10] . They assume 
that is an effective hard-sphere radius and predict 
that the value of the rms displacement at the fee transi- 
tion, Atrans, Satisfies 

2(1 - t:V2/6)-^/^ Atrans - Rws - (14) 

where Rws — (47r/3)~^/^a is the Wigner-Seitz radius. 
They then use an approximate formula for the effective 
frequency in an Einstein-like model to determine A for 
the fee and bcc structures. These values of A give two 
predictions for the transition line. As shown in Figure |S1 
their predictions are qualitatively correct but lie roughly 
10-40% above 

The discrepancy in Fig. |3 could be due to a failure ei- 
ther of Eq. (|14|l or of the approximations used to find 
A. To test this we performed equilibrium simulations 
at ft{.;*„^ in both bcc {N = 3456) and fee {N = 4000) 
systems. Our results for the rms displacements, A fee 
and Abcc, are compared to the predictions of Eq. ifTH) in 
Table IVTll The rms displacements for the fee structure 
lie quite close to the prediction for small A, and about 
13% above it at A = 7. Finite size effects decrease A fee 
relative to the N = cx3 value [5| . These results indi- 
cate that most of the error in Vaulina et. al.'s transition 
lines comes from substantial underestimation of the rms 
displacements. They calculate A in the harmonic ap- 
proximation, and anharmonic corrections increase A for 
these A [IJ. Note that the measured bcc displacements 
in Table VI are larger due to the bcc lattice's softer shear 
modes [2l|. 

V. SUMMARY AND CONCLUSIONS 

We calculated the bcc-fcc coexistence line of Yukawa 
systems to an uncertainty of approximately 1% through 
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TABLE VII: Rms displacements at T{^1„^{X). Atrans is the 
prediction from Ea ll4l while A fee and Atcc are results from 
equilibrium MD simulations. 

X Atrans Afcc A^cc 



3 0.0915 0.089 0.096 

4 0.1181 0.118 0.128 

5 0.1341 0.140 0.154 

6 0.1447 0.159 0.177 

7 0.1522 0.171 0.192 



integration of the mechanical work along Bain transfor- 
mation paths. The range of bee stability was found to 
be slightly greater than that found in previous compre- 
hensive studies [^[12,123, and the triple point shown to 
lie at higher inverse screening length. The large changes 
in Tfrans witli Tc for small A indicate that the relative 
stability of fee and bcc phases depends sensitively on 
long-range correlations, and calls into question the use of 
local nearest-neighbor arguments to calculate the tran- 
sition line. Nevertheless, we found that one such phe- 
nomenological criterion jlO|, derived from the idea that 
the fee phase is stable when interparticle interactions are 
hard-sphere- like j predicts the transition line remark- 
ably well when combined with MD results for the mean- 
squared displacement. 

Comparison with lattice-dynamics results shows that 
anharmonic terms in the total energy favor the bcc phase 
for all A, but that these corrections are overbalanced by 
anharmonic contributions to the entropy. The change in 
entropy appears to reflect an increase in the frequency of 
long-wavelength shear modes in the bcc phase. This in- 
crease also stabilizes the bcc phase against a linear shear 



instability observed for A > 7.67 at low T. 

We found that shifts in the transition line due to finite 
size effects are less than 1% if N-3000-4000, but that 
the presence of long-range order at temperatures near the 
transition line in weakly screened systems requires a cut- 
off radius larger than that used in some previous studies 
Accurate simulations of weakly screened (A < 3) 
systems in this temperature range require either larger 
system sizes or an Ewald-like summation over periodic 
images 23]. However, we have also shown that a reason- 
ably small potential cutoff need not introduce large errors 
in a transition line calculation in the moderate-screening 
regime, provided the cutoff is chosen with some care. 

It is known that the phase behavior of real systems 
such as charge-stabilized colloidal suspensions is not fully 
described by pointlike Yukawa interactions. Recent sim- 
ulations of charged macroions in a dynamic neutralizing 
background have shown that the repulsive interactions 
between macroions are truncated by many-body effects, 
destabilizing the crystalline phases in the weak screening 
limit li^ 03' Independently, hard-core repulsions 
significantly alter the phase diagram when the volume 
fraction is more than a few percent [48ll4^l50l |. However, 
we hope that our high-precision calculation of the point 
Yukawa fcc-bcc transition line may serve as a benchmark 
for further studies of more sophisticated models. 
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